Multi-component lattice-Boltzmann model with interparticle interaction 

Xiao wen Shan 12 and Gary Doolen 2 

A previously proposed [X. Shan and H. Chen, Phys. Rev. E 47, 1815, (1993)] lattice 
Boltzmann model for simulating fluids with multiple components and interparticle forces 
is described in detail. Macroscopic equations governing the motion of each component are 
derived by using Chapman-Enskog method. The mutual diffusivity in a binary mixture 
is calculated analytically and confirmed by numerical simulation. The diffusivity is 
generally a function of the concentrations of the two components but independent of the 
fluid velocity so that the diffusion is Galilean invariant. The analytically calculated shear 
kinematic viscosity of this model is also confirmed numerically. 
KEY WORDS: lattice-Boltzmann; multi-phase flow; diffusion. 

1. Introduction 

Recently the lattice-Boltzmann Equation (LBE) method has emerged as a new promising method of 
computational fluid mechanics (CFD). This method was developed from a discretized fluid model 
known as the Lattice Gas Automaton (LGA).( 1>2 ) In the LGA model, fluid is modeled microscopically 
as a collection of particles moving on a regular lattice along the links. The particles collide with each 
other on lattice sites according to some carefully designed collision rules which conserve the number 
of particles and momentum. The coarse-grained fluid variables, such as density and fluid velocity, 
can be shown to obey a set of macroscopic equations which are very similar to the Navier-Stokes 
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equations. Since only a small number of bits are required to characterize the states of each lattice site, 
and the collision operation is local in most cases, the LGA can be implemented very efficiently for a 
large number of lattice sites. Fluid motion can then be simulated without integrating partial differential 
equations. 

Although the LGA method has a few highly desirable advantages over the conventional methods, 
it also suffers from intrinsic drawbacks. For instance, the large statistical fluctuation limits the practical 
usage of the LGA method. To suppress the statistical noise, McNamara and Zanettit 3 ) suggested 
modeling the LGA with a lattice Boltzmann equation. Mean population, instead of the discrete 
particles, are used to simulate fluid flows. A similar approach was proposed by Higuera et al( 4 ) in 
which a linearized collision term is used in the lattice Boltzmann equation so that the algorithm is 
computationally more efficient and can be easily generalized to three dimensions. More recently, it 
was pointed out( 5 6 ) that in addition to the suppression of statistical noise, the unphysical artifacts in the 
original LGA, known as the lack of Galilean invariance and the velocity-dependent pressure term, can 
also be eliminated when a single relaxation time collision term (also known as the BGK collision term 
after Bhatnagar, Gross and Krook( 7 )) with a proper choice of the equilibrium distribution function is 
adopted in the lattice Boltzmann equation. Comparisons with conventional CFD methods demonstrated 
that such a LBE model can give quite satisfactory results in simulating both hydrodynamics and 
magnetohydrodynamics problems.^ 8 ' 9 ) 

An interesting and important application of the LGA/LBE methods is the simulation of fluid flows 
with interfaces between multiple phases. There are numerous complex flow systems in both natural and 
industrial processes that involve convection-coupled mass transfer near fluid interfaces. The density 
gradients in such problems are often so large that the conventional linear diffusion equation cease to be 
an acceptable approximation. Such problems have posed considerable difficulties to the conventional 
CFD methods, especially when the interfaces can undergo topological changes. Since the formation 
of fluid interfaces is microscopically due to the long-range interaction between the molecules of the 
fluid/ 10 ) the separation and reconnection of the interfaces require additional terms to be inserted in 
the Navier-Stokes equations. While the method of Molecular Dynamics can treat interfacial problems 
by taking into account the details of the inter-molecular interaction, it becomes computationally too 
expansive when the large scale flow structures have to be modeled simultaneously. 
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Because the LGA and LBE methods were developed based on microscopic description of fluids, 
interactions between fluid elements can be naturally included. Several LGA/LBE models for multi- 
phase flows have been developed since the first introduction of the LGA. Rothman and Keller^ 11 ) 
developed the first LGA model for simulation of two immiscible fluids. A Boltzmann equation version 
was formulated later.t 12 ) The same idea was used to achieve a low diffusivity in the simulation of 
partially miscible fluids^ 13 ) with LBE method. In these models, a repulsive force between the two 
fluid components is most intense in the interfacial zone where the "color gradient" is large. Since the 
"interfacial zone" is defined by an arbitrary small parameter, for the LBE model, a finite-amplitude 
perturbation is required for a phase separation to occur in a initially homogeneous system. The early 
stages of phase separation can be expected to have some arbitrariness. In addition, these models 
have other problems such as the anisotropy of the surface tension^ 14,12 ) and difficulties in dealing with 
components with different densities. Appert et al( 15 ) suggested another LGA model to simulate a 
liquid-gas type of phase transition. Attractive or repulsive forces exist between particles that are several 
lattice units apart. As the range of the force increases, the system separates into two phases. However, 
the low efficiency and all the problems that plague the LGA methods have prevented this model from 
being used to solve practical problems. 

In a previous paper,( 16 ) we proposed a scheme to incorporate the interparticle forces in LBE 
models with multiple components. Interaction potentials of different nature are defined between 
particles of the same or different components. Each component has its own molecular mass and 
relaxation time. The LBE model was shown to be able to simulate a fluid with an arbitrary non-ideal 
gas equation of state. When the equation of state is properly chosen, phase separation occurs both in 
single- and multiple-component systems. In most cases, the interaction can be restricted to involve only 
nearest neighbors so that the model is computationally efficient. Isotropy of the surface tension and 
the density profile across the liquid-gas interface have been shown both analytically and numerically in 
a single component fluid containing two phases. ( 17 ) In order for this model to be used in quantitative 
simulations of multiphase flow problems, the densities of the two phases have also been obtained 
analytically as functions of a temperature-like parameter. 

In a system with more than one component, there could be a very complicated dependence on 
the details of the interaction potential. Several phases with different densities and concentrations and 
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multiple diffusion processes exist in such a system. In general, the gradients of the concentrations 
can be very large and the diffusivity depends on the concentration of each component. Before the 
LBE model can be used to simulate multiphase flow with mass transfer between phases, the diffusion 
process has to be well understood and the diffusivity calculated. In this article, we give the details of the 
multiple component LBE model with interparticle forces. In Section 2, we derive the macroscopic fluid 
equations governing the motion of each component using the Chapman-Enskog expansion method. We 
obtain the mutual diffusivity in a binary mixture as a function of the concentrations of the components. 
It does not depend on the fluid velocity. In Section 3, the diffusivity in a numerical experiment on 
a two-dimensional hexagonal lattice was measured and found to be in excellent agreement with the 
analytical values. Problems involving convection-diffusion processes near a two-liquid interface and 
evaporation in porous media have been investigated using the present LBE model. These simulation 
results will be discussed in future publications. 



2. Multiple component Lattice Boltzmann model 

We now review the multiple-component LBE model( 16 ) with interparticle forces. Consider the motion 
of particles of S different components on a regular lattice in D -dimensional space. The particles of 
the <7th component have the molecular mass of m CT , a = 1, . . . , S. The population of the particles of 
the <7th component having the velocity e a at lattice site x and time t is denoted by tj£(x, t); where, 
{e a ; a = 1, . . . , b} is the set of vectors of length c pointing from x to its b neighboring sites. The 
evolution of ti£(x, t) is described by the following lattice Boltzmann equations: 

<(x + e a ,f+l)-<(x,0=--k(x,0-< (eq) (x,Ol; a=l,---,S. (1) 

where r CT is the collision interval of the <7th component. On the right hand side, we adopt the BGK 
collision term and chose the equilibrium distribution functions to be n£ ea ^ = N a (n a , u^), where 



N a (n, u) = < 



(1^*0 ,D_ , D(D+2) . _ Du^_\ . — i ... K 

U \h+ r2h e °- U + 2c 4 fc eaea • UU 2c 2 b J ' a ~ A > >° n . 



n[do-K); 



4 



do > is an arbitrary constant. It was shown that this choice of the equilibrium distribution corrects 
the unphysical artifacts of the LGA and yields the correct Navier-Stokes equations. ( 5 ) In the equations 
above, = J2 a n a * s tne num ber density of the <rth component, are parameters which will be 
calculated from the distribution functions and the long-range interactions. 

In the absence of interparticle forces, all the components are ideal gases. We assume that the 
velocities u^, a = 1, . . . , S, all equal a common velocity u'. Since in this case the total momentum of 
particles of all components should be conserved by the collision operator at each lattice site, it follows 
directly from Eq. (1) that 

u' = £—/£-. (3) 

where p a = m^n^ is the <7th component mass density, and p^u^ = Y^ a n Z e a is the <rth component 
momentum. 

We let the strength of the long range force between particles of component a at site x and 
particles of component a at site x' be proportional to the product of their "effective masses," 
V> (T (ti ( 7-(x))V'o-(tIo-(x')). The "effective mass" of the <7th component ^(^(x)) is defined as a function 
of the local density n a . The form of ip^n) can be arbitrarily chosen and will determine the equations of 
state of the fluid components and composite fluid.t 17 ) Summing over all the components and interacting 
sites, the total long range force acting on the particles of the <7th component at site x is: 

F CT (x) = - W(x) J2 E G -( x > x')^(x')(x' - x), (4) 

x' o- 

where the Green's function satisfies G^x, x') = G^x', x). If only homogeneous isotropic interac- 
tions between the nearest neighbors are allowed, G^x, x') reduces to the following symmetric matrix 
with constant elements, 

0, |x - x'| > c 



G^(x,x') = < 



(5) 

Q a „, |x - x'| < C. 



Furthermore, if the densities vary at a spatial scale much larger than the lattice spacing, c, F CT can 
be approximated by: 

We limit the discussion in this paper to interactions which are homogeneous, isotropic, and between 



nearest neighbors. 

The long-range interparticle force introduced above causes an extra momentum change to the 
<rth component in addition to the momentum exchange caused by collisions with other components. To 
incorporate this momentum change in the dynamics of the distribution functions, one can simply define 

p^u? = / o CT u' + r CT F CT . (7) 

It can be verified that at each time step, the momentum change of the <7th component due to long range 
interaction is F CT . 

When a long-range force exists, the collision operator does not locally conserve the total mo- 
mentum of all the components, although it does conserve the number of particles of each component. 
By summing Eq. (1) over all directions, we obtain a mass-conservation relation for each of the S com- 
ponents. Multiplying Eq. (1) by e a and m CT , and summing over all the b directions and S components, 
we can obtain the change of total momentum at each site. These relations are: 

£<(x + e a ,*+l)-£<(x,*) = 0; (7=1,..., S, (8) 

a a 

J2 m a J2 <( x + e a, * + l)e a - E m " E <( x > *) e ° = E F - < 9 ) 

a a a a a 

The second equation states that during a collision, the total momentum of the particles at a lattice site 
is changed by the interactions with the particles on neighboring sites. Nevertheless, the momentum of 
the whole system can be shown to be conserved. We find the change of total momentum of the whole 
system by summing over all the lattice sites, 

AP = E F -( x ) = " E E 9™ E Vv(x)^(x + e a )e a . (10) 

X <7 X a,a a 

Here, e a is a dummy variable. If there is no net momentum flux at the boundaries, as in the case with 
periodic boundaries, x can be viewed as a dummy variable too. After changes of variable e a — > -e a 
and then x - e a — > x, we have 

A p = EE^EW x )^( x + e ^ e - ( u ) 

X (T,& CL 
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Since the matrix Q ad . is symmetric, we have AP = -AP, and therefore AP = 0. 

The question arises as how to calculate the macroscopic "fluid velocity" from the distribution 
functions since the momentum of the <rth component are m CT J2 a n°e a before a collision and (1 - 
Y^ a n a e a + 7^ S a n a^e a after collision, as obtained by summing the distribution function 
before and after collision. They differ by a significant amount in region containing large density 
gradients such as the interface between two phases. Since it is the average of the two values that 
represents the over-all mass transfer, we define the velocity of the whole fluid, u, to be given by the 
average of the two momentum values summed over all the components. We can show that this fluid 
velocity vanishes at equilibrium. Using Eqs. (7) and (3), we find the average to be 

/0U = ^m CT ^<e a + ^F CT , (12) 

a a <7 

where p = p a is the total mass density of the fluid. The second term is generally negligible except 
in the interfacial zone. 

We now follow the Chapman-Enskog method of successive approximation^ 1819 ) to obtain the 
macroscopic fluid equation of the multiple-component LBE model. The equilibrium distribution 
functions, n*, are expanded as an infinite series, and the time-derivative d/dtis also divided into parts 
accordingly: 

r=0 r=l 

Since now the collision operator depends upon the spatial derivatives of the densities, it is impossible to 
choose a distribution function which is function of velocity and density only and makes the collision term 
vanish. We therefore let the leading order distribution functions of all components be the equilibrium 
distribution about the fluid velocity, u, namely n^ ^ = N a (n a ,u),so that a set of macroscopic equations 
in terms of the correct fluid variables can be obtained as the result of the expansion procedure. The 
Boltzmann equation will be satisfied at the next order when terms that depend on spatial derivatives are 
included. Because £) a n^ ^ = n a and "V n a^e a = pu, the next order terms will satisfy the 
following equations: 

£< (1) = o; E^E< (1) ^ = -^E F - (W) 



Taylor-expanding the left-hand side of Eq. (1) to second order and using the expansions (13), we 



have 



<(x + e OJ f+l)-<(x,f) 



a °-(°) a °-(°) a <r( [ ) 

di n a ' , _ 0.(0) , d 2 n a K ' , d[Ti a ' 
- e a ■ vn a 



dt 



dt 



dt 



On substituting Eq. (15) into the conservation relations Eqs. (8)-(9) and collecting terms of leading 
order, the following equations are obtained: 



ding. 
~dT 
dipu 
dt 



V • (n a u) = 0, a 
~c\\ - d ) 



D 



p + /OUU 



l,--,S 



(16) 
(17) 



Using Eq. (16), and noting that Y^a^o- — -fs^Ew ^W*> Eq. (17) can be written in the 
following form 

D,n 1 

(18) 



Dt 



■-Vp, 
P 



where D[/Dt = di/dt + u • V, and p is the pressure given as function of the densities of all the 
components by the following equation of state 



(l - d )p + - Ga*i>*i>* 



(19) 



The second second term in Eq. (19), which causes the equation of state of the fluid to be one of a 
non-ideal gas, depends explicitly on the interparticle force. 

The terms of the next order in the expansion of the mass-conservation relation yield the following 
equation for each component: 



dip* 
dt 



m CT V • £ <We a + ^V.(|E < (0) e a + V • £ <(°)e a e a ) = 0. (20) 

a \ a a / 
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Using Eq. (16), we can evaluate the third term. 



^f|E< (0) ^ + v -E< (0)e ^) 

\ a a / 



Pv C 2 (l - dp) 

— Vp + V/v- 

P D 



(21) 



To calculate J2 a m o- n a^e a in the second term of (21), we substitute the expansions of Eq (13) up to 
second order into the kinetic equation (1), obtaining the following equations: 



Jft E m o-< (0) e a + V • ^2 m a n^e a e a = 



^(u - u') + ^m CT < (1) e a 



F CT , (22) 



which when multiplied by and summed over all the components, become 



P (u - u') = £t ctFct + \ y: f ct + y, - e^ ct . 



(23) 



In the derivation above, Eqs. (14) and (21) are used. The expression £) a m a ria^e a can then be solved 
from Eq. (22) in terms of the macroscopic fluid variables. Combining Eq. (20) with the first order 
equation (16), we obtain the following equations at second order 



^ + V • (p„u) = -r CT V • F CT + (t„ - l -\ V 



c 2 (l-d ) 

= Vpo- 

D p 



E t.F. + ^ E ^ + \ E - E ^ 



(24) 



When summed over all the components, the equations above become the continuity equation of the 
whole fluid at second order: 



dp 
dt 



V • (pu) = 0. 



(25) 



Mutual diffusivity in a two-component system can be calculated using Eq. (24). For notational 
convenience, we rescale the interaction constants Q a5 . as G a „ = bQ a ^l{\ - do). Using Eq. (19), 
Eq. (24) can be written as 



dpi 
dt 



V-Goiu) 



c\l-d {) ) 



D 



c\\-d Q ) 



D 



V • {AVpi - BV P2 ) 

V • (BV P2 - AW pi) . 



(26) 
(27) 
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Here A and B are functions of the densities of the two components 



. pifpiTi + pm 1\ , ., p\T 2 + P2T\ G\ \1pip2 - Gntpipi ,„ ON 

A = — +^1 ( 2 8) 

P V P 2 / /> /> 

„ PlfPlT2 + P2Tl 1\ , ., p\T 2 + p 2 TlG 22 ^ 2 p\- G\ 2 ^\p 2 ,- n , 

5 = — +^2 , (29) 



where = and ^ = dip^ / dp a \ 

We consider the evolution of small perturbations of the densities. Let p a = p® + p\ + ■ ■ where 
p® are uniform, time-independent equilibrium densities of the two components, and p^ <C p® are small 
perturbations. Define the concentration of the <rth component as c CT = p®/p°, where p° = p\ + p\. 
Obviously c\ + C2 = 1. We linearize Eqs. (26)-(27) about the equilibrium densities to obtain 

^ + .?V-u = ^^}v • (A°V P \ - B«V P \) (30) 
^ + ^V-u = ^^V.(^-^V,I) (31) 

where D/Dt = d/dt + u • V, and 

A° = C 2 (ciT 2 + C 2 Ti - i) + lf)[(ciT 2 + C 2 n)(G n lf)lC 2 - Gl2^2Cl) (32) 
5° = Ci(ciT 2 + C 2 Ti - i) + V4(C1T2 + C 2 Ti)(G 2 2V'2Cl - G\2^\C 2 ). (33) 

In a pure diffusion process, if we assume the sound speed is much larger than the speed of 
diffusion, the pressure field can be taken as uniform. When the equation of state is not an ideal gas one, 
the total density field is not a uniform constant during the diffusion. The density perturbations p\ are 
therefore not independent of each other, and a velocity field with a generally non-zero divergence will 
be generated. In the present case, p\ and p\ are related by Eq. (19): 

[1 + i>[{G n i>i + G 2 vh)} P \+[l + ^{Gn^i + G 22 i>2)] pi = 0. (34) 

We can then write p\ and p\ in terms of a single perturbation field, £(x), 

/>J(x) = [1 + ^(GWi + G 22 $ 2 )\ £(x) (35) 
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= -[l + il>i{Gwl>i + G2nh)]Z(x)- (36) 
It follows from the continuity equation (25) that 

[^{GnTpi + G221P2) - i>[{Gni>i + G 2 iih)] ^ + p°V ■ u = 0. (37) 

On substituting the three equations above into either Eq. (30) or (31), we find that £(x) satisfies the 
following diffusion equation 



D£ = c 2 (l - dp) 
Dt D 



1 

(C1T 2 + C 2 Ti)7 - - 



V 2 £, (38) 



where, the effect of long range interaction on the diffusivity is absorbed into the single factor 

= (1 + G n M[)(l + G 22 ^ 2 ) ~ gj^MljM4 (39) 

1 + ClGnV'lV'l + C2G22^2^' 2 + G 12 (ciV'l^ + c 2^2^l) ' 

7 becomes unity in the absence of long range interaction. Furthermore, if the relaxation times of the 
two components are equal, namely t\ = r 2 = r, the diffusivity becomes c 2 (l - c?o)(r - 1/2)/ D, in 
agreement with other authors/ 13 ' 20 ) 

A few remarks are called for at this point upon the diffusivity derived above. First, the diffusivity 
depends on the relative density concentration of the two components as expected. This dependence, 
through not only the c CT 's but also the ipo- 's, is quite general and can theoretically be tuned to approximate 
any given diffusivity. We are therefore able to simulate a class of diffusion problems in which the 
diffusivity can vary and depend on concentrations. Second, unlike the LBE model for miscible fluid 
due to Flekk0y,( 20 ) in which the mutual diffusivity depends on the velocity of the fluid, the diffusivity 
does not have a velocity dependence. The diffusion in the present LBE model is fully Galilean invariant. 

When applied to convection-diffusion problems, we notice that the diffusivity can be varied inde- 
pendently of the shear viscosity by changing the long range interaction potential. Generally attractive 
forces between particles of same component, and repulsive forces between different components both 
tend to decrease the diffusivity. 

The second order momentum equation can be derived in a way similar to that of one-component 
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fluid without interparticle interaction^ 19 ) It turns out that long range interaction does not affect the 
shear viscosity, at least to second order. When multiple components are present, the shear viscosity is 
simply v = c 2 (£ CT c CT r CT - ±)/(Z> + 2). 

3. Numerical simulation 

In this section, we present results of numerical simulations of the LBE model described above. We first 
measure the mutual diffusivity in a binary mixture by studying the decay of a sinusoidal concentration 
wave with small amplitude. The simulation was carried out on a 256 x 16 hexagonal lattice. Initially 
the perturbations of the densities of the two components were set up according to Eq. (34) so that the 
pressure was uniform in the whole field. The relaxation times of the two components are t\ = 0.7 
and T2 = 1.0. Both number densities were set equal to unity. The "effective masses" were chosen 
to be if>i(ni) = 1 - exp(— Tii) and if>2 = ^2- The diffusivity is obtained by measuring the decay rate 
of the concentration wave. The results of the measurement are plotted in Fig. 1 as function of Gn, 
which measures the strength of the force between particles of the two components. G\\ and G22 are 
fixed at positive numbers 0.01 and 0.02 respectively, so that interactions between all the particles are 
repulsive. Two sets of data are presented for two different values of the molecular mass of the second 
component, mi = 2 and mi = 4. The first component always has a molecular mass of 1.0. The 
density concentration ci is then 0.33 and 0.2 in the two cases. In Fig. 1, the solid lines are the analytic 
predictions, which are in excellent agreement with the experimental values. 

The mutual diffusivity is also measured in the presence of a uniform velocity field parallel to the 
concentration gradient. No dependence of the diffusivity on the velocity is observed. For instance, 
when all the parameters are chosen to be the same as in the case above in which 7712 = 2 and G12 = 0.04, 
the diffusivity measured without flow is 0.0549, compared with 0.0550 when a uniform flow of speed 
0.2 is present. 

The expression for the shear viscosity is also verified similarly by measuring the decay rate of a 
sinusoidal transverse wave on the same lattice. The concentrations of the two components are varied 
by changing the molecular mass ratio while keeping the total density of the fluid constant. Without 
loss of generality, we chose the number densities 711 = ni = 1 and p = 5 in the simulation. The 
relaxation times of the two components are 0.6 and 1.2 respectively. Gu equals 0.04 here and all other 
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parameters are the same as before. The measured viscosity is plotted in Fig. 2 as function of the density 
concentration of the first component. It is again in excellent agreement with the analytic prediction 
indicated by the solid straight line. 

4. Conclusions 

We have described in detail a multiple-component LBE model with interparticle interactions. The 
equations governing the evolution of the densities of the components are derived using the Chapman- 
Enskog technique. The mutual diffusivity and shear viscosity in a binary mixture are derived analytically 
using those equations and the agreement is obtained to high accuracy with numerical simulations. The 
diffusivity generally has an explicit dependence on the concentrations of the two components and can 
be tuned to model physical diffusivities which are concentration-dependent. When the interparticle 
interaction is turned off and the relaxation times of the two components are equal, the diffusivity does 
not depend on the concentrations, in accordance with other authors. The Galilean invariance of the 
diffusivity in this model is also confirmed by numerical simulations. 

In the current paper, we focus our attention to the derivation of macroscopic equations and the 
transport coefficients of multiple-component system, which can be fully or partially miscible. In the 
latter case, this model can be used to simulate fluid flows which involve phase transition and interfaces 
between different phases. The derivation of the transport coefficients in this case becomes more valuable 
because in some cases the concentration gradients near an interface are so large that the diffusivity can 
not be approximated as a constant. The dependence of the diffusivity on the concentration itself 
is essential in the study of convection-diffusion problems that involve interface. In principal, the 
phase-diagram, the density-profile across an interface, and the surface tension can be calculated for a 
multiple-component system in a way similar to, but more complex than, a one-component system.^ 17 ) 
We defer the detailed discussion to a future publication. 
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Figures 



Fig. 1 . Mutual diffusivity in a binary mixture, as function of the strength of the repulsive force between 
the two components. Solid lines are theoretical predictions and the symbols are the results of 
numerical simulation. The two sets data are for two different density ratios of the two components. 

Fig. 2. Shear viscosity as function of the density of the first component. Interparticle interactions exist 
in the simulation but this has no effect on the shear viscosity. The solid line is the theoretical 
prediction. 
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